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Parametric study of the interface behavior between two immiscible 
liquids flowing through a porous medium 

Alejandro David MariottiP* Elena Brandaleze, 2 *" Gustavo C. BuscagliaP 

When two immiscible liquids that coexist inside a porous medium are drained through 
an opening, a complex flow takes place in which the interface between the liquids moves, 
tilts and bends. The interface profiles depend on the physical properties of the liquids 
and on the velocity at which they are extracted. If the drainage flow rate, the liquids 
volume fraction in the drainage flow and the physical properties of the liquids are known, 
the interface angle in the immediate vicinity of the outlet (0) can be determined. In this 
work, we define four nondimensional parameters that rule the fluid dynamical problem 
and, by means of a numerical parametric analysis, an equation to predict is developed. 
The equation is verified through several numerical assessments in which the parameters 
are modified simultaneously and arbitrarily. In addition, the qualitative influence of each 
nondimensional parameter on the interface shape is reported. 



I. Introduction 

The fluid dynamics of the flow of two immiscible 
liquids through a porous medium plays a key role 
in several engineering processes. Usually, though 
the interest is focused on the extraction of one of 
the liquids, the simultaneous extraction of both liq- 
uids is necessary. This is the case of oil production 
and of ironmaking. The water injection method 
used in oil production consists of injecting water 
back into the reservoir, usually to increase pressure 



* E-mail: mariotti.david@gmail.com 

1" E-mail: ebrandaleze@frsn.utn.edu.ar 

* E-mail: gustavo.buscaglia@icmc.usp.br 

1 Instituto Balseiro, 8400 San Carlos de Bariloche, Ar- 
gentina. 

2 Departamento de Metalurgia, Universidad Tecnologica 
Nacional Facultad Regional San Nicolas, 2900 San 
Nicolas, Argentina. 

3 Instituto de Ciencias Matematicas e de Computacao, 
Universidade de Sao Paulo, 13560-970 Sao Carlos, Brasil. 



and thereby stimulate production. Normally, just 
a small percentage of the oil in a reservoir can be 
extracted, but water injection increases that per- 
centage and maintains the production rate of the 
reservoir over a longer period of time. The water 
displaces the oil from the reservoir and pushes it to- 
wards an oil production well [1]. In the steel indus- 
try, this multiphase phenomenon occurs inside the 
blast furnace hearth, in which the porous medium 
consists of coke particles. The slag and pig iron 
are stratified in the hearth and, periodically, they 
are drained through a lateral orifice. The under- 
standing of this flow is crucial for the proper design 
and management of the blast furnace hearth [2] . In 
both examples above, when the liquids are drained, 
a complex flow takes place in which the interface 
between the liquids moves, tilts and bends. 

Numerical simulation of multiphase flows in 
porous media is focused mainly in upscaling meth- 
ods, aimed at solving for large scale features of in- 
terest in such a way as to model the effect of the 
small scale features [3-5]. Other authors [6-8] use 
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the numerical methods to model the complex mul- 
tiphase flow that takes place at the pore scale. 

In this work, we numerically study the macro- 
scopic behavior of the interface between two im- 
miscible liquids flowing through a porous medium 
when they are drained through an opening. The 
effect of gravity on this phenomenon is considered. 
We define four nondimensional parameters that 
rule the fluid dynamical problem and, by means 
of a numerical parametric analysis, an equation to 
predict the interface tilt in the vicinity of the orifice 
(0) is developed. The equation is verified through 
several numerical cases where the parameters are 
varied simultaneously and arbitrarily. In addition, 
the qualitative influence of each non-dimensional 
parameter on the interface shape is reported. 

II. Parametric Study 

The numerical studies in this work were carried out 
by means of the program FLUENT 6.3.26. Differ- 
ent models to simulate the two-dimensional para- 
metric study were used. The volume of fluid (VOF) 
method was chosen to treat the interface problem 
[9]. The drag force in the porous medium was 
modeled by means of the source term suggested by 
Forchheimer [10]. The source term for the i th mo- 
mentum equation is: 



= - + \pC\f\Vt 
\ a 2 



(i) 



For the constants a and C in Eq. (1), we use the 
values proposed by Ergun [10]: 



3 d 2 
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C = 1.75 
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(2) 



(3) 



where e is the porosity, d is the particle equivalent 
diameter, p is the density, V is the velocity and p 
is the dynamic molecular viscosity. 

Considering that the subscript 1 and 2 represent 
the fluid 1 and the fluid 2 respectively, three nondi- 
mensional parameters were considered in the para- 
metric study: viscosity ratio, \ir = pi/p2, density 
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Figure 1: Sketch of the numerical 2D domain. 

ratio, pR = p\j p<i-> and nondimensional velocity, 
Vr = V0P2L/P2] where Vb is the outlet velocity 
and L a reference length. 

i. Domain description 

The numerical domain considered to carry out the 
parametric study was a two-dimensional one com- 
posed by the porous medium sub-domain and the 
outlet sub-domain. The porous sub-domain is a 
rectangle 10m wide and 10m tall. Inside of it, a 
rigid, isotropic and homogeneous porous medium 
was arranged. We use a porosity and particle di- 
ameter of 0.32 and 0.006m, respectively. 

For the outlet domain we use a rectangle 0.02m 
wide and with a height L = 0.01m divided into 
two equal parts and located at the center of one 
of the lateral edges. The part located at the end 
of the outlet domain is used to impose the outlet 
velocity. Figure 1 shows a complete description of 
the domain. 

A quadrilateral mesh with 2.2 x 10 4 cells was 
used, where the outlet sub-domain mesh consists 
of 200 elements in all the cases studied. 

As boundary conditions, on edge 1 we define a 
zero gauge pressure condition normal to the bound- 
ary and impose that only the fluid 1 can enter to the 
domain through it. On edge 2 the boundary condi- 
tion is the same as on edge 1 but the fluid consider 
in this case is fluid 2. On edge 7 we impose a zero 
gauge pressure normal to the boundary but in this 
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Figure 2: Interface evolution when the initial posi- 
tion is below the outlet level, without gravity. 

case the fluids can only leave the domain. On the 
other edges (edges 4, 5, 6, and 8) we impose a wall 
condition where the normal and tangential velocity 
is zero except for edge 3, at which the tangential 
velocity is free and the stress tangential to the edge 
is zero. 

ii. Interface evolution 

To illustrate how an interface reaches the stationary 
position from an initially horizontal one, three sets 
of curves were obtained. 

Figure 2 shows the interface evolution for the 
case without the gravity effect and the interface 
initial position is below the outlet level. The inter- 
face modifies its tilt to reach the exit and it changes 
its shape to reach the stationary profile. 

Figures 3 and 4 show the interface evolution 
when the gravity is present but the interface initial 
position is below and above the outlet level respec- 
tively. 

iii. Viscosity effect 

One of the most important parameters to modify is 
the dynamical viscosity of fluid 1. We maintain the 
properties of the fluid 2 as the properties of water 
(density 998Kg/m 3 , and dynamical viscosity 0.001 
Pa.s) and the density of fluid 1 as the density of the 
oil (850 Kg/m 3 ). The dynamical viscosity of fluid 
1 was varied from values smaller than those of fluid 
2, to values much greater. 
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Figure 3: Interface evolution when the initial posi- 
tion is below the outlet, with gravity. 
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Figure 4: Interface evolution when the initial posi- 
tion is above the outlet, with gravity. 

Two sets of curves were obtained, one consid- 
ering the effect of gravity and the other without 
considering it. Figure 5 shows the stationary in- 
terface profiles for the different values of viscosity, 
without gravity. A value of Vr = 1.5 x 10 4 and 
Pr — 1.17 were chosen. It is possible to observe 
that, when fluid 1 has a viscosity higher than that 
of fluid 2, the interface profile is above the outlet 
and points downwards at the outlet. If fluid 1 has 
a lower viscosity, the opposite happens. 

When considering gravity, the value of Vr was 
changed to 1 x 10 5 (Vo = lOm/s), since for smaller 
values the interface may not reach the outlet (this is 
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Figure 6: Stationary interface profiles modifying 
the fluid 1 viscosity with the gravity effect. 

later studied in Fig. 10). Figure 6 shows the curves 
obtained for this situation, where the interface only 
lies over the outlet level for the higher {ir values. 

iv. Outlet velocity effect 

Vo is varied from a small value, similar to the porous 
medium velocity (Vo = 0.2m/s or Vr — 2000), to a 
very large one (Vo = 50m/s or Vr = 5x105). Main- 
taining the properties of fluid 2 similar to those of 
water, two sets of curves were obtained (fiR > 1 
and \i r < 1), shown in Figs. 7 and 8, respectively. 

When the effect of gravity was considered, two 
additional sets of curves (Figs. 9 and 10) were ob- 



Figure 8: Stationary interface profiles for several 
values of Vr, without gravity, for [ir < 1 (fiR = 
0.01). 

tained. 

Figure 7 shows the effect of Vr when the viscosity 
of fluid 1 is greater than that of fluid 2, without 
gravity. It is possible to see that, as Vr increases, 
the interface tilt at the outlet is maximal for Vr = 
1 x 10 5 . 

On the other hand, Fig. 8 shows the interface 
profiles when the viscosity of fluid 1 is smaller than 
that of fluid 2. We observe that as Vr increases the 
interface tends to the horizontal position. 

Figure 9 shows the different stationary interface 
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Figure 9: Stationary interface profiles for several 
values of Vr, with gravity, for pr > 1 (pr = 35). 
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Figure 11: Stationary interface profiles for several 
values of the density of fluid 1, with Vr = 1.5 4 and 
p R 35. 



fluid 2 with the properties of water and the viscos- 
ity of fluid 1 as 0.035Pa.s (pr = 35), the density 
of fluid 1 was varied from its original value to one 
three times smaller than that of fluid 2. In Fig. 11 
it is seen that as pR increases, the interface profile 
ascends significantly, with a less significant change 
in the tilt angle at the outlet. 



III. Generic expression 



Figure 10: Stationary interface profiles for several 
values of Vr, with gravity, for pr < 1 (pr = 0.01). 

positions when the gravity effect is present for pr > 
1. The effect of gravity is quite significant, the 
interface ascends but only for the highest value of 
Vr it lies above the outlet level. 

Figure 10 shows the curves when the viscosity of 
fluid 1 is lower than that of fluid 2 (pr < 1). The 
behavior is different from that without gravity. In 
fact, there exists a minimum outlet velocity below 
which the interface does not reach the outlet. 

v. Density effect 

The effect of the density ratio on the interface pro- 
file was studied in the presence of gravity. Keeping 



From the study on the influence of each nondi- 
mensional parameter on the interface behavior, an 
equation that predicts the interface angle at the im- 
mediate vicinity of the outlet (9) was crafted. For 
practical reasons, the cases where gravity is present 
were considered to develop the equation. 

In Sect. 2.2, it is possible to see that when the 
nondimensional parameters pR, pr and Vr are con- 
stant the interface changes its shape until it reaches 
a stationary profile. 

For this reason, a fourth nondimensional param- 
eter is considered, the volume fraction of fluid 1 in 
the outlet flow (VF). 

A generic expression [(Eq. (4)], consisting of 
three terms and containing 22 constants, was ad- 
justed by trial and error until satisfactory agree- 
ment with the numerical results was found. 
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3.2 
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-1.26 











Table 1: Constant values in the generic expression. 



PT 


D 


£ 


1/a 


C 


Resistance 


A 


0.005 


0.3 


10.88 x 10 Y 


1.81 x 10 4 


Very high 


B 


0.006 


0.32 


5.88 x 10 7 


1.21 x 10 4 


High 


C 


0.02 


0.2 


3 x 10 7 


1.75 x 10 4 


Medium 


D 


0.02 


0.25 


1.35 x 10 7 


8400 


Low 


E 


0.05 


0.17 


8.41 x 10 6 


1.18 x 10 4 


Very low 



Case 


Fluid 1 


Fluid 2 


Mi 


M2 


Pi 


P2 


1 


Heavy 
oil 


Water 
solu- 
tion 


0.4 


0.005 


850 


998 


2 


Light 
oil 


Water 
emul- 
sion 


0.012 


0.06 


850 


998 


3 






0.08 


0.008 


4000 


4680 


4 


sene 


Water 


24 x 10 -4 


0.001 


780 


998 


5 


Ace- 
tone 


Water 


3.3 x 10~ 4 


0.001 


791 


998 


6 






0.003 


0.01 


400 


720 


7 






0.2 


0.01 


1500 


2700 


8 






0.3 


0.004 


3300 


5940 


9 


Light 
slag 


Hot 

Pig 
iron 


0.02 


0.001 


2800 


7000 


10 






0.5 


0.013 


500 


1250 


11 


Medium 
slag 


Pig 
iron 


0.4 


0.005 


2800 


7000 


12-16 






0.08 


0.008 


4000 


4680 


17-21 


Heavy 
slag 


Pig 
iron 


0.4 


0.005 


2800 


7000 



Table 2: Porous medium types. 



= a&V&i 
+ WrKPr eM-i^p l R VF m ) 
+ nu° R V& R eM-r^ R Vy R VFn (4) 

Table 1 shows the values of the constants in the 
generic expression. 

i. Equation verification 

To verify that the generic expression (4) predicts 
the value of correctly when the parameters are 
arbitrarily modified, 21 additional numerical cases 
were simulated. These cases cover a wide range of 
physical properties of the liquids and of the char- 
acteristics of the porous medium (by means of the 
coefficients 1/a and C). 

The interface angle obtained from the generic 
expression (Oqe) was compared with the interface 
angle obtained from the simulations (Os). Table 2 
shows the five porous medium types (PT) that were 
chosen. 

Some cases (1, 2, 4, 5, 9, 11, 17-21) were chosen 
based on the possible combinations of immiscible 
liquids that can be manipulated in real situations. 
For the remaining cases, the properties of the liq- 
uids were fixed at arbitrary values (fictitious liq- 
uids) so that a wide range of the nondimensional 



Table 3: Cases description. 



Case 


PT 


PR 


PR 


V R 


VF 


1 


B 


1.17 


80 


5 x 10 4 


2.8 


2 


B 


1.17 


0.2 


5 x 10 4 


87.2 


3 


B 


1.17 


10 


10 x 10 4 


28.4 


4 


B 


1.28 


2.4 


5 x 10 4 


96.1 


5 


B 


1.26 


0.33 


1 x 10 5 


96.8 


6 


B 


1.8 


0.3 


1 x 10 5 


93.2 


7 


B 


1.8 


20 


1 x 10 5 


16.1 


8 


B 


1.8 


80 


8 x 10 4 


18.4 


9 


B 


2.5 


20 


1 x 10 5 


100.0 


10 


B 


2.5 


40 


1 x 10 5 


5.5 


11 


B 


2.5 


80 


5 x 10 4 


70.8 


12 


A 


1.2 


10.0 


1 x 10 5 


23.1 


13 


B 


1.2 


10.0 


1 x 10 5 


28.4 


14 


C 


1.2 


10.0 


1 x 10 5 


41.5 


15 


D 


1.2 


10.0 


1 x 10 5 


53.4 


16 


E 


1.2 


10.0 


1 x 10 5 


63.7 


17 


A 


2.5 


80.0 


5 x 10 4 


15.8 


18 


B 


2.5 


80.0 


5 x 10 4 


24.3 


19 


C 


2.5 


80.0 


5 x 10 4 


43.2 


20 


D 


2.5 


80.0 


5 x 10 4 


70.8 


21 


E 


2.5 


80.0 


5 x 10 4 


94.0 



Table 4: Parameter values and PT for all cases 
described in Table 3. 

parameters was covered. Table 3 shows the descrip- 
tion of each numerical case, while Table 4 shows 
the corresponding nondimensional parameter val- 
ues and PT. 
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Figure 12: Comparison between the predictions of 
the generic expression and the numerical result for 
the 21 validation cases. 



We define an error (e = 1OO|A0|/18O) as the 
percentage of the absolute value of the difference 
between the interface angles (A6 = 0s — @ge) di- 
vided by the interface angle range (180°). Figure 
12 shows the comparison between the generic ex- 
pression and the numerical cases. It is seen that 
the generic expression (4) predicts the interface an- 
gle, for the cases used in this study, with an error 
smaller than 10%. 



IV. Conclusions 

A numerical study of the macroscopic interface 
behavior between two immiscible liquids flowing 
through a porous medium, when they are drained 
through an opening, has been reported. Four 
nondimensional parameters that rule the fluid- 
dynamical problem were identified. Thereby, a nu- 
merical parametric analysis was developed where 
the qualitative observation of the resulting inter- 
face profiles contributes to the understanding of the 
effect of each parameter. In addition, a generic ex- 
pression to predict the interface angle in the imme- 
diate vicinity of the outlet opening (0) was devel- 
oped. To verify that the generic equation predicts 
the value of correctly, 21 numerical cases with 
widely different parameters were simulated. Con- 
sidering that the cases encompass a large class of 
liquids and porous media, the prediction of within 
an error of 10% is considered satisfactory. 
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